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It is well-known that molecular dynamics integrators, which are used for lattice quantum chro- 
modynamics (QCD), suffer from instabilities and possess a rather low order of the accuracy. 
Hence, it is highly desirable to construct a new class of geometric integrators, that overcomes 
these instability problems and increases the order of accuracy without increasing remarkably the 
computational costs. 

Sin this paper we consider for this purpose multistep methods and give an overview of known 
, results to systematize important knowledge for such methods being the right choice for lattice 

QCD simulations. At the end we try to answer the question: can multistep method be used as 
molecular dynamic integrators and what might be the advantage of it. 
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1. Introduction 

In this paper we will give a short introduction to symplectic multistep integrators, review some 
recent results and comment on the possible later application in lattice QCD computations. 

A geometric integrator is a numerical method that preserves the geometric properties of the 
exact flow of an autonomous ordinary differential equation (ODE) 

y=/(y), y^K^. (1.1) 

Especially, when solving numerically a Hamiltonian problem, it is of paramount importance that 
the chosen scheme retains some properties of the underlying continuous problem, like the time 
reversal symmetry or the area preservation property. 

We consider a compatible linear multistep methods (LMM) with k steps of the form 

k k 

^ apn+j ^hY, Pjf{yn+j)^ (1.2) 

where h = At denotes the time step of the grid tj = + jK and aj, Pj are real parameters, aj ^ 0, 
and |ofo| + |j8o| > 0. For an application of (1.2) we need an initial value yo = y(^o) as well as starting 
approximations , . . . ^yk-i to 3;(/o) , • • • ^y{tk-i), that are usually obtained by Runge-Kutta methods. 
We recall that the LMM (1.2) is of order s if and only if (cf. [8]) 

t^; = 0, ta/ = ^ti8/-\ l<i<s, t^^,/+V(^+l)ti8,/. (1.3) 

J=0 j=0 j=0 j=0 j=0 

Furthermore, the multistep method (1.3) has the two characteristic polynomials 

p(^) = ta,^^ cT(^) = ti8,-^^ (1.4) 

and the LMM (1.2) is called irreducible if these polynomials (1.4) have no conmion roots. 

Since the stability analysis for LMMs was difficult to establish Dahlquist [2] suggested to 
consider instead so-called one-leg methods that only need one evaluation of the forcing function /. 
The one-leg method (OLM) associated to the multistep method (1.2) is defined by 




(1.5) 



where we assume the normalization condition 0"(1) = Y!j=oPj = L The LMMs (1.2) and OEMs 
(1.5) are closely related. Let 3; denote the sequence of approximated values obtained from a LMM 
(1.2) and y be the sequence of approximated values obtained from a OLM (1.5), then we have 

k 

Hence the analysis of the stability for LMMs can be reduced to the stability analysis of OEMs. 
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In the sequel we focus on the case that the ODE (1.1) is a linear autonomous Hamiltonian 
system, i.e. p = 2n (where n is the number of degrees of freedom in mechanics) and 

y = J-WH{y), y G where / = ( ° | . (1.6) 

Here, H : denotes a smooth Hamiltonian function. Now we can rewrite the LMM (1.2) 

k k 

£ ajyn+j ^hf, ^jJ-'VH{yn+j). (1.7) 



2. Area-Preservation and Time-Reversibility of Multistep Methods 

A numerical integrator is used in a molecular dynamics step of the Hybrid Monte Carlo al- 
gorithm has to satisfy an area-preservation property, which follows from the symplecticity of the 
numerical method, and a time-reversibility property, which is an extension of the symmetry. 

The symmetry is fulfilled if the coefficients of a scheme (1.2) (or (1.5)) satisfy the relations 

ak-j = -aj, Pk-j = Pj for all j = 0, 1, . . . 

i.e. 

They are also charaterized by having an odd number of time steps, i.e. ^ = 2v — 1, v=l,2, 

This means that the numerical solutions satisfy the following reversibility requirement: whenever 
y^, . . . ,yn+k satisfy the relation (1.2) (or (1.5)), yn+k^ ---^yn satisfy (1.2) with h replaced by —h. 
From this it follows that symmetric multistep methods are time-reversible. 

Definition 1. A mapping g : M?^ M?^ is called symplectic (with respect to J) if 



\dg{y)] 


T 

J 


\dg{y)] 


[ dy \ 




[ dy \ 



It is well-known, that the solution of the ODE (1.6) at any fixed time regarded as a function 
on the initial data 3;(0) G M?^, (so-called phase flow) is a symplectic mapping. Hence, it is a natural 
task to seek for numerical methods that retain this property (in a sense to be specified later). 

In the literature there exist (at least) two definitions for symplectic LMMs. Eirola and Sanz- 
Serna [9] considered the transformation g : R^"^^ M^"^^, Y£ F^+i, where Y£ = {yj, . . . 
that is associated to the LMM (1.7), and obtained the following positive result. 

Theorem 1 ([9]). Assume that the OLM is symmetric and irreducible. Then the corresponding 
mapping ^ F^+i is symplectic with respect to the matrix J, where A — {Xij) is given by 

^ij — ^ {(^i+mPj+m + ^^y+mft+m)? i ^ 0, J < k. 

m>0 
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But here the symplecticity is regarded with respect to another skew symmetric matrix A /. 
We note that this symplecticity is equivalent to the preservation of quadratic first integrals, cf. [1]. 
For example, Ge and Feng [6] showed that the standard second order leap-frog scheme 



jn+i-yn^^hJ ^WH{yn+i) is symplcctic wxt. ( with A 





Contrary, in a different second approach. Tang [10] proved a negative result for the step- 
transition operator (underlying one-step method) introduced by Feng [5] g : M?^ satisfying 

k k 

t^ajg^^hf^^jfog^, (2.2) 

7=0 j=0 

where g^ stands for j-time composition g'.gogO"'Og. This operator characterizes the LMM (1.2) 
as y\ = g(yo), • --^yk^ §{yk-i) = <?^(yo), • • • , e.g. for Hamiltonian systems the LMM (1.7) reads 

k k 

Y^ajg^^hY^^jJ-\WH)ogK (2.3) 

This step-transition operator g allows for a definition of symplecticity for LMMs: 

Definition 2 ([10]). The LMM (1.7) is symplectic if the operator g defined by (2.3) is symplectic, 
i.e. (2.1) holds for all Hamiltonian functions H and all sufficiently small step-sizes h. 

Theorem 2 ((Conjecture of Feng) [10]). For a Hamiltonian system (1.6), any compatible LMM 
(with order s>\ of form (1.7) is not symplectic. 

Despite this negative result of Theorem 2, it is known that the second order mid-point rule 

yn+\ -yn = hJ-^VH{]^{yn+i (2.4) 
is a symplectic multistep method. This fact motivates to consider generalized LMMs of the form 

k k k k 

t,ajyj^hf^^jJ-'VH{f^yjiyi), with f^Jji^l^ J = 0,...,L (2.5) 

j=0 j=0 /=0 /=0 

Unfortunately, also for the scheme (2.5) the result is rather negative. 

Theorem 3 ([10]). For a Hamiltonian system (1.6), if a difference scheme (with order s>\ of the 
form (2.5) is symplectic, then it must be of order 2. 

In fact, Dai and Tang [3] showed that (2.4) is the only symplectic scheme of this form. How- 
ever, following the concept of G-stability proposed by Dahlquist there exists another third way of 
transfering the definition of symplecticity to multistep methods, cf. [2]. 

However, recall that for integrators to be suitable for molecular dynamics integration, it is 
sufficient to satisfy the following area-preservation property. 

Definition 3. A mapping g : M?^ E^" is called area-preserving (without an orientation) if 

. dg{y) 



det 



dy 

which is a shghtly weaker assumption than simplecticity 



= 1. (2.6) 
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3. Numerical Experiments 



We consider a model of a simple harmonic oscillator (SHO) with the Hamiltonian 

H{p,q) = ^p^ + ^w\^ 



(3.1) 



to investigate the stability behavior of solutions for this system by G-symplectic multistep methods. 




^'^^^* **** *** y ** 
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Figure 1: Solutions of the SHO; explicit Euler method with step size /z = 0.1, initial value 
{po^ ^o) = (0, 1); implicit Euler method with step size /z = 0.1, initial value (po, ^o) = (0, 1). 

We choose the three G-symplectic methods and show the corresponding numerical solution of 
the SHO problem (3.1). The first method is a 4-step explicit method 



h 



(3.2) 



The second method is a predictor-corrector method 



Jn+A - yn+?> ^ (55/^+3 - 59/^+2 + 37/^+ 1 - 9fn) as a predictor 
h 

Jn+A - yn+?> ^ (9A+4 + 19/^+3 fn+l) as a corrcctor. 



(3.3) 



using the explicit method to compute /^+4 of the implicit one. Finally, we consider a partitioned 
method, where for each equation of the system of SHO we apply a different multistep method: 



yn+?> - yn+2 + yn+l —yn^h {fn+2 + fn+l] 
yn+3 - yn+l = 2/z {fn+2 + fn+l ) • 



(3.4) 



Figures 2, 3, 4 show (from left to right) the numerical solutions obtained by the corresponding 
multistep methods (3.2) - (3.4) and their long-time behavior {n = 1000000 time steps). 
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Figure 3: Solution of the SHO by the predictor-corrector method (3.3) with step size /z = 0. 1, initial 
value (po, ^o) = (0, 1); the error of the computed solution of the method (3.3) to exact solution. 




Figure 4: Solution of the SHO by the partitioned method (3.4) with step size /z = 0.1, initial value 
{po^ ^o) = (0, 1); the error of the computed solution of the method (3.4) to exact solution. 

The method (3.2) shows quite stable behavior for the numerical solution of the SHO, but 
further considerations show that increasing the number of steps lead us to oscillations and sym- 
plecticity of the solutions will be destroyed. But for the short term problems the method behavior 
is suitable. The method (3.3) collapses after some time and the error of this method grows expo- 
nentially, but still for short time period it gives the proper results and the symplectic property is 
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satisfied. The last method (3.4) yields the best results, its numerical solutions behave symplec- 
tic even during long time, it conserves the energy of system properly and this class of partitioned 
multistep method, according to [8], gives the correct results even for the long time integrations. 

4. Conclusion and Outlook 

In spite of the collected negative results, the numerical experiments showed that the short-time 
behavior of the multistep methods is rather promising. The main advantage of multistep method is 
that the high-order version of such methods can be easily obtained by one function evaluation per 
time step and it, consequently, will increase the accuracy of computations. 

Despite these predominant negative statements, we recall that in lattice QCD simulations one 
solely needs an area-preserving integrator which is a slightly weaker assumption than the discussed 
symplecticity. In a forthcoming paper we will investigate, following an idea of Hairer [7], a pro- 
jected LMM that conserves the Hamiltonian and hence makes the acceptance step in the hybrid 
Monte Carlo simulations obsolete. This feature is especially interesting for small lattice spacings. 
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